Fast Simulation of Facilitated Spin Models 
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We show how to apply the absorbing Markov chain Monte Carlo algorithm of Novotny to simulate 
kinetically constrained models of glasses. We consider in detail one-spin facilitated models, such as 
the East model and its generalizations to arbitrary dimensions. We investigate how to maximise the 
efficiency of the algorithms, and show that simulation times can be improved on standard continuous 
time Monte Carlo by several orders of magnitude. We illustrate the method with equilibrium and 
aging results. These include a study of relaxation times in the East model for dimensions d = 1 to 
d — 13, which provides further evidence that the hierarchical relaxation in this model is present in 
all dimensions. We discuss how the method can be applied to other kinetically constrained models. 



I. INTRODUCTION 

By nature, glassy systems are dynamically slow, with 
relaxation times found to increase rapidly as the temper- 
ature T is lowered Q, 0, • O ne consequence of this slow 
down is the difficulty in probing the dynamics of such sys- 
tems in the long time regime through means of numerical 
simulation. Traditional numerical methods often become 
insufficient to obtain data within a realistic time frame 
thanks to time scales growing at best exponentially with 
inverse temperature and in most cases much faster. A 
common algorithm used for systems close to dynamical 
arrest is the rejection free algorithm known as the "n- 
fold way", or continuous time (CT) Monte Carlo 0,13 . 
While CT can provide an impressive improvement over 
standard Monte Carlo (MC), it can still become ineffi- 
cient for some extremely slow systems. 

The natural generalization of the CT algorithm is 
Novotny's Monte Carlo with Absorbing Markov Chains 
(MCAMC) 6], which so far has been mainly used to 
study magnetic reversal. In this paper we show how to 
apply the MCAMC method to simulate kinetically con- 
strained models (KCM) of glass formers Q, in particular, 
facilitated spin models such as the Fredrickson- Andersen 
(FA) model [8| and the East model Q in arbitrary di- 
mensions. We show that in the case of East models 
the MCAMC method allows to improve simulation times, 
both for equilibrium and out-of-equilibrium dynamics, by 
several orders of magnitude at low temperatures. 

The paper is organized as follows. In Section [n] we 
outline the MCAMC technique. Section 11111 provides a 
detailed analysis of its application to the one-dimensional 
East model. In Section ITVl we discuss several important 
approximations which allow to maximise the computa- 
tional gain. The method is generalized to the East model 
in arbitrary dimensions in Section^] and, in Section IVll 
to models that interpolate between the East and FA cases 
[Toj . Section lYlII describes faster, higher level versions of 
MCAMC. In Section lVlIII we present speed tests of the al- 
gorithms. In Section llXl we show some illustrative results 
obtained with the MCAMC method, including a study of 
the T dependence of relaxation times in the East model 
for dimensions d = 1 to d = 13. We conclude in Section 
13 discussing the possible implementation of the method 



to simulate other KCMs. 



II. OUTLINE OF THE MCAMC METHOD 

Continuous time Monte Carlo works well for systems 
where the vast majority of moves are most likely to be re- 
jected due to kinetic or energy constraints In order to 
avoid constantly attempting and failing to make a move, 
as in standard MC, in CT one fast-forwards the system, 
always accepting moves, and updating the clock by how 
long it would have taken in standard MC. This can speed 
things up greatly in systems with slow dynamics because 
making an unfavourable move can take a long time. 

The trouble with most slow systems is that often after 
an unlikely move, even in the CT scheme, the most likely 
move is to undo what has just been done before. With the 
MCAMC method of Novotny @, instead of just fast for- 
warding to when one is accepted, one can fast-forward to 
when two (or more) unlikely moves have been accepted, 
updating the clock appropriately. To do this whilst keep- 
ing the dynamics exact and keeping to detailed balance 
one must use the formalism of absorbing Markov chains 
@ (see ^3 f° r a pedagogical review) . 

In an MC algorithm any given move depends only on 
the two states that the system is moving between and not 
on any previous moves. This is by definition a Markov 
process and allows us to treat the MC algorithm as a 
Markov chain 12]. A Markov chain is characterised by 
the matrix M which defines transition probabilities be- 
tween states. If the vector x T (m) indicates the proba- 
bility distribution of the system after iteration to, the 
probability distribution at the next step to + 1 is given 
by x 1 '(to + 1) = 2? t (to)M. 

We define an absorbing Markov chain by separating 
the available states into s transient states and r absorb- 
ing ones 0, ^J. The system always starts in a tran- 
sient state and by successive applications of the Markov 
matrix explores the transient subspace until it lands in 
an absorbing (or exit) state from where it cannot leave. 
We can divide the general state vector x T into absorbing 
and transient parts, to get the (r + s)-dimensional vector 
x T = (u r ,v r ) where contains the transient states. 
The initial state in this form must obey xf = (0 T ,vf). 
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With this structure the Markov matrix can be written in 
the form, 



III. APPLICATION OF MCAMC TO 
KINETIC ALLY CONSTRAINED MODELS 



m = ( i rxr r xs 

"■s X r J- s X s j 



(1) 



where I is the identity matrix, is the zero matrix and 
subscripts indicate the size of each matrix. The posi- 
tions of the identity and zero matrices guarantee that if 
the system falls into an absorbing state then it does not 
leave. The transient matrix, T, gives the probabilities 
for moving between transient states and the recursive 
matrix, R, gives the probabilities for moving from the 
transient states to the absorbing states. 

For a given starting vector vf the probability of still 
being in the transient subspace, ptrans., after m steps is 



Pit 



(2) 



where e is a vector with all elements equal to 1. The prob- 
ability of absorbing to a particular state after m steps is 
given by summing over the probabilities of absorbing at 
each time step. This gives the vector p[ hs , 



Pahs, after m 



I + T- 



)R. 



(3) 



If the exit has taken place at step m, then the proba- 
bilities of absorbing into the different exit states is given 
by: 



-T 

Pubs .at m 
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(4) 



Here it is convenient to introduce the fundamental ma- 
trix 



N = (I-Tp 1 =I + T + T 2 



(5) 



which can be used to obtain the probability that the sys- 
tem will absorb to a particular state irrespective of when 
it exits, 



Plbs. = NR. 



(0) 



The fundamental matrix can also be used to determine 
the average time to leave the transient subspace 0, 



(r) = vjNe. 



(7) 



Once our system is in the initial state we can generate 
an exit time by solving the inequality 



Vj 1 e < r < Vj 1 e, 



(8) 



where r is a random number between and 1. Next, 
we use a second random number to choose an absorp- 
tion state from the distribution in equation 0} and then 
we can update the system appropriately. The new state 
will become the initial state in another absorbing Markov 
chain, and so on. A successful MCAMC algorithm will 
choose transient states such that the system tends to 
move between them many times before exiting. 



The MCAMC algorithm has been used in the study 
of magnetic reversal and nucleation in systems such as 
the Ising model, see e.g. [1,0,0]. In the case of mag- 
netic reversal, for example, there is a well defined initial 
transient state corresponding to the metastable config- 
uration in which all spins are aligned in one direction. 
This state represents a deep minimum in the energy: on 
reversing a single isolated spin the overwhelming likeli- 
hood is that the system will immediately relax back to 
the initial state. An MCAMC algorithm with two tran- 
sient states (s — 2) overcomes this problem: the system 
can escape from the metastable configuration by insist- 
ing that two spins are flipped simultaneously; this process 
is characterised by two transient states, the metastable 
configuration and that corresponding to a single spin flip 

6]- 

Kinetically constrained models (KCMs) are frequently 
used to study the dynamical behaviour of glass form- 
ers 0- Of particular interest is the analysis of relax- 
ation times and the emergence of length-scales . Sim- 
ilarly to the magnetization reversal problem above, at low 
temperatures or high densities, these systems evolve by 
falling into deep energy or free energy traps from which 
it is difficult to escape. The difference, however, is that 
there is a multiplicity of trapping states, not just one as in 
the Ising model case, which change as the system evolves 
dynamically. In order to apply the MCAMC method to 
KCMs one needs to identify the nature of the transient 
states in these systems. 

As an example consider the East model in one- 
dimension 0, 13 • The East model consists of a chain 
of Ising spins, rii — {0; 1}, upon which a directional facil- 
itation constraint is imposed: a given spin rii may only 
flip if its nearest neighbour to the left is in the excited 
state, rii— i = 1- The rates for allowed moves are such 
that the equilibrium distribution is that corresponding 
to the energy function H = J J^. rii (we set J = I from 
now on), 
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11, 11 



10, 



(9) 



where e = e _/3 and j3 = 1/T. The result is a model 
in which excitations propagate in the eastward direction. 
The East model has been found to reproduce the dynamic 
behaviour of fragile glassy materials 0, • 

From the rate equations above one can see that it is 
energetically unfavourable to excite a down spin the pro- 
cess becoming increasingly unlikely as temperature is de- 
creased. Consider a configuration with only isolated ex- 
citations: 

•■■100---100---100--- 

When in this state the only choice is to excite a spin and 
pay the accompanying energy penalty. Consequently it is 
highly likely that the raised spin is immediately relaxed, 
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hence returning to the previous state. This is analogous 
to the all-up configuration of the Ising model example. It 
is important to note that in this isolated state all excita- 
tions must be separated by a minimum of two unexcited 
spins. When separated by only a single spin there are 
two possible outcomes from the creation of an excita- 
tion, i.e. one can create a double state 110 or a triplet 
111. This produces an unnecessary complication since 
the algorithm can no longer be classified by two simple 
transient states. 

By analogy to the Ising model one may define two tran- 
sient states for the system, the isolated configuration de- 
scribed above and states in which a single excitation pair 
exists. However, unlike the case of magnetic reversal, it 
is clear that neither of the transient states identified for 
the East model are unique. It is possible to construct 
numerous configurations which satisfy the above crite- 
rion, in essence we have identified two classes of tran- 
sient state. Once again the absorbing states consist of all 
configurations attainable by the excitation of two spins, 
either forming two isolated doubles or a triplet state, 

• ■■110- •■110 •■■100- • • 

•■■111---100---100--- 

For the East model it is possible to classify each lattice 
site according to its local neigbourhood. Taking a site 
along with its nearest and next-nearest neighbour to the 
right, each site can be classed according to a binary la- 
belling scheme, i.e. 100 = 4, 110 = 5, etc., where the 
number of sites in each class is, etc. Using this 

notation we define the entry condition for the algorithm 
with s = 2 transient states as the point at which the 
number of sites in class 4 equals the total number of ex- 
citations present in the lattice, M, i.e. N4 = M. 

Before constructing the transient and recursive ma- 
trices it is necessary to determine the probabilities for 
all possible transitions between the different states. The 
transient and recursive states may be labelled as follows, 



100- 


• 100 • • • 


Vl 


110- 


■ 100 ■ ■ ■ 


V 2 
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• 100 • • • 


a 2 



with the following transition probabilities 
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where N is the system size. 

These transition probabilities are then used to build 
the transient and recursive matrices for the system 



T 
R 



1 — x x 
V 1-x-y 



x - ey ey 



(10) 
(11) 



where x 



— ana y — . 

The absorption probabilities for the ui and u 2 states 
can be found by taking the fundamental matrix, N, and 
solving equation © giving 



P( Ul ) 

P{U2) 



1 - 
1 



1 



(12) 
(13) 



where we have used an initial state vector iff = (1 0). 

To determine the exit time one must choose a random 
number and iteratively solve the inequality given in equa- 
tion ijHJ). One then proceeds to choose an exit state from 
the distribution formed by the exit probabilities, equa- 
tions (|12fl and l|13|) . It is clear that matrices T and R 
are characterised by the variable 7V4 and as such both 
the probability distribution for the absorption states and 
the exit time are governed by the entry state, each state 
having its own unique solution. 

This s = 2 construction provides an algorithm that 
improves on standard continuous time, s = 1, by a fac- 
tor proportional to This improvement in com- 
putational speed is offset by the algorithmic complexity 
required to formulate the s = 2 model. 



IV. APPROXIMATIONS FOR THE UPDATE 
TIME 

Computationally, the most expensive part of the algo- 
rithm as described above is the procedure used to deter- 
mine the time to exit from the transient state. To per- 
form the calculation exactly involves diagonalising the 
T matrix and iteratively solving the inequality using the 
halving method 5] or something similar. There are, how- 
ever, a number of approximations that we can employ to 
get around this. The exact form for vfT m e for the s = 2 
case is 



A™ + XT - (X^ - XT) ( ^ + z 



where 



+ 4eN4 and Ai, A 2 



. (14) 
are the eigenval- 



ues of T. Both eigenvalues are quite close to (and less 
than) 1. However, in the limit of large m, we have that 
(A2/A1)'" <C 1, allowing us to simplify equation (fT4*jl . If 
we drop the restriction that m must be discrete then we 
can write © as an equality, 



log 



2r 

1 + z + 1/4? 



log(Ai), 



(15) 
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where again r is a random number between and 1. Both 
z and Ai depend on N4 and can be stored in a lookup 
table. The approximation works best when m is large, 
so for low temperatures (T < 1) where the time steps 
are larger it works very well. At higher temperatures one 
must be careful using this approach. 

Another possibility is to free oneself from the require- 
ment to pick the update time from a distribution and use 
instead the average. This does mark a departure from the 
exact Monte Carlo algorithm, but in most cases it turns 
out to be a reasonable simplification (it is analogous to 
the approximation made when going from the n-fold al- 
gorithm Q to the CT one 5]). If we take the average 
time, then we can use equation JJJ which requires calcu- 
lation of the fundamental matrix, N, either analytically 
or numerically. For the East model, N only depends on 
the number of excitations M, so the time updates can be 
stored in a lookup table allowing for a significant increase 
in speed. 

To check the validity of using the average value for time 
updates instead of picking them from a distribution, we 
can use the result 

(r 2 ) = vj (2N 2 -N)e, (16) 

with J7| to calculate the mean square fluctuations. This 
shows that for lower temperatures the error on any given 
measurement is ~ (r). Whilst this seems large it is im- 
portant to remember that we are always looking at log- 
arithmic time and on this axis the error is less signifi- 
cant. Also there are many iterations between sampling 
points and the measurements are averaged over many 
runs which will help to reduce any discrepancy. All the 
simulations for this paper were performed using the av- 
erage time update. 
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xd — ey eyd 



(18) 



It is easy to solve for the exit probabilities to each of the 
absorbing states, where A4 again indicates the number 
of excitations in the isolated state, 



P(ui) = 1- 
P{u 2 ) = 



d - 1 + N 4 d 
d 



d-l + N 4 d 



(19) 
(20) 



The average lifetime to exit from the transient subspace 
may be obtained from l(7jl. 



e 2/3 + (2N A d + d- 1)^ 

N 4 d(N 4 d + d-l) 



(21) 



For systems in which (3 is high and N4 is large com- 
pared to d, one may approximate the exit time to 



=2/3 



(N 4 d) 2 ' 



(22) 



hence time steps are a factor of d 2 smaller than those of 
the East model in d = 1. 



VI. FA-EAST CROSSOVER MODEL 

The MCAMC algorithm described above can also be 
applied to the FA model 0, 0] and, more generally, to 
the model that interpolates between the FA and East 
models 0] 1 which serves as a simple model for fragile- 
to-strong transitions. This model is characterised by the 
rates 



V. GENERALISATION TO ANY DIMENSION 

The method described in the previous section can eas- 
ily be extended allowing one to construct generalised 
transient and recursive matrices for the East model in 
any spatial dimension d. Considering the transient states 
for the system it is clear that the s = 2 algorithm is trig- 
gered when all excitations within the lattice are isolated 
by a region of space which encompasses all moves attain- 
able by two successive spin flips. The d = 2 analog of the 
"100" above is: 



100 







1 



where no triangles may overlap if the algorithm is to trig- 
ger. As for the d = 1 case, the T and R matrices are ob- 
tained by evaluating the probabilities for all transitions 
between the transient and absorbing states. This analysis 
yields matrices of the following form 



1 — xd 
V 



1 — y — xd — (d — l)ey 



(17) 
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01, 01 



be 



11, 



11 ^ 10, 



10 



(l-6)e 



11. 



(23) 

The limit b — > corresponds to the East model, and 
the limit b — > 1/2, to the FA model. For intermediate 
values of b the model displays a crossover between East- 
like dynamics at higher temperature, to FA dynamics at 
low temperature. 

The MCAMC algorithm is applied in much the same 
way as in the East model case, except that now we have 
to allow for the possibility of movement to the west. In 
the simplest version, the transient states are the same 
as for the East model, and the absorption states are in- 
creased to include any move to the left. The result is an 
s — 2, r — 5 absorbing Markov chain, with the following 
transient states 



0100- ■ -0100 • 



Vl, 



0110- • -0100- ■• v 2 , 
and absorbing states 

0110---0110--- mi, 



5 



0111 •••0100- •• 



«2, 



transient state. In one dimension one can represent the 
states as follows, 



1100 ■■■0100 ■ 



100- • -100- ■ -100- 



vi, 



0110- •• 1100- •• or 1110- ■ -0100- •■ u 4 , 



■110---100---100--' 



>'2- 



0100 ■••0100 • 



Us. 



The two states in m are the same as far as this algorithm 
is concerned. Caution would be required if they were to 
be used as transient states. The transient and recursive 
matrices are as follows 



1 — xd xda 
ya 1 — y — xd — a(d — l)ey 



R 







bxd 



a(xd — ey) aedy bxd by 



(24) 



(25) 



where a = 1 — b, and x and y are the same as for the East 
model case. 

From here on the procedure is exactly the same as 
with the East model except that there are more absorb- 
ing states to choose from. One may use Eq. JHJl to ob- 
tain values for the absorption probabilities. Caution is 
required as the approximation breaks down in the regime 
of high temperature and high symmetry (6 — > 1/2). It 
should be noted that, while less striking, even in the FA 
limit of b = 1/2 this algorithm outperforms standard CT. 



VII. HIGHER ORDER MCAMC 

The entirely isolated state is problematic in terms of 
the dynamics of the East model. In order to relax the 
isolated excitations must propagate in the lattice until 
they encounter another excitation along the direction of 
facilitation. Movement of this nature is promoted by the 
occurrence of branching events, 



100 



111 -> 101. 



The "triplet" absorption state, U2, is the rate limiting 
step for branching events, and hence the propagation of 
excitations in the lattice. However, from the absorption 
probabilities, Eqs. I|12l) and (|13[) . we see that compared 
to the u\ state the u 2 exit state is suppressed by a factor 
of N 1 _ 1 - Hence, the formation of triplets is unlikely, 
particularly when the system size is large. 

To overcome this problem it is possible to extend the 
MCAMC algorithm to include the u\ state as a transient 
state of the system. There are now three transient and 
three absorbing states, one absorbing state corresponds 
to the u 2 state of the s = 2 algorithm, the remainder 
corresponding to configurations attainable from the new 



110- •■110- • -100- 



110- • -110- • -no- 



111 •••100- • -100- ■ 



• 111 ••■110- • -100- • 



V3, 



iti, 



U 2 , 



lt 3 . 



Once again the transient and recursive matrices may 
be constructed by considering all possible transitions be- 
tween the states. 
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1 — x x 
y 1 — x — y x — ey 
2y l-2y- 


ey 
x — 2ey 2ey 



(26) 
(27) 



The s = 2 transient matrix, Eq. (TDJ, is now a subma- 
trix of T; the addition of an extra transient state has 
appended one extra row and column to the matrix, the 
rest of the structure remaining intact. 

Unlike the case of s = 2, it is not so simple to generalise 
the s = 3 matrices for any dimension. This arises from 
the non-equivalence of the i>3 state in dimensions d > 1, 
i.e. in two dimensions 



10 / 1 •■■10 
110 10 10 

When considered as absorption states the two configu- 
rations above may be treated identically since the prob- 
ability of exiting to each state is the same. However, 
as transient states each configuration has different exit 
probabilities and as such must be treated independently. 
In essence, one requires an s = 4 algorithm to provide 
the equivalent result in dimension two and above. 

Returning to the d = 1 example, we find that the u 2 is 
now the most likely absorption state. This is because all 
other exit states require the excitation of an additional 
spin, i.e., they are suppressed by a factor of e° . Solving 
for the average lifetime gives 

p 2/3 
M W — . 

v 1 N 4 

Hence, s = 3 improves on s = 2 by a factor of N4. While 
this may seem a modest enhancement in performance, 
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note that the extra algorithmic complexity required to 
develop s — 3 is negligible. Having made a working s = 2 
algorithm one may essentially use s = 3 for free. 

As CT enables one to obtain an e' 3 speed increase over 
traditional MC, s — 2 enables one to achieve a further tu 
improvement of e' 3 over CT. In effect, s — 2 enables one .5 
to bypass all e' 3 processes (i.e. those that involve the 
excitation of a single spin) by insisting that two succes- 
sive spins are excited. The double and triplet states of 
the s = 2 model are examples of e 2 ^ 3 processes. In or- 
der to construct an algorithm with a further e' 3 speed 
gain requires one to identify all of the e 2 ' 3 arrangements 
and include them as transient states. This means that 
the absorption states now correspond to all configura- 
tions attainable from the transient states which result in 
the simultaneous excitation of three spins. This analysis 
leads to an s = 7 algorithm consisting of seven transient 
and seven absorbing states. 

In one dimension, s = 7 may be triggered when all ex- 
citations within the lattice are separated by at least three 
unexcited spins, i.e. 1000 • • • 1000. To maximise perfor- 
mance it is useful to use a hybrid algorithm consisting 
of s = 1, 3 and 7 components with each sub-algorithm 
activated by its own triggering condition. 

In order to improve algorithmic efficiency it is conve- jq^ 
nicnt to compute absorption probabilities using Eq. © 
rather than the exact form of Eq. (0}. Unlike the case 
of s = 2, where it may be shown that two expressions 
are identical, for higher order algorithms the solution of 
Eq. © only provides an approximation. In general one 
must employ caution when using this approach. For both 
s = 3 and s = 7 it has been shown that the approxi- 
mation is good for all regimes in which the algorithms 
are effective, the approximation breaking down at higher 10 
temperatures. 



VIII. SPEED TESTS 
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FIG. 1: Temperature dependence of CPU time for equi- 
librium East model trajectories of total Monte Carlo time 
t = 10 7 x e 2f3 and system size N = 10 5 , for MC, CT, and 
MCAMC algorithms. The straight lines indicate the approxi- 
mate speed-up of the MCAMC simulations. CPU time shown 
relative to the average time needed when using CT dynamics. 
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In the low T limit, the average exit time from an 
s = 2 MCAMC algorithm iteration for the East model 
is approximately e 2 ° /(Nid) 2 . The corresponding aver- 
age time step for standard CT is e@ /N±d. The s = 2 
time step becomes larger by a factor of e@ /N±d. It gets a 
speed-up from e^ 3 , and a slowdown from as the more 
excitations that are present upon entering the algorithm 
the quicker it exits, and from d, as the higher the di- 
mension the more facilitated sites are available. At low 
temperature, however, JV4 « Ne, so for fixed system size 
N the speed-up factor of s = 2 MCAMC with respect to 
CT grows as e 2/3 . 

In figures 1 and 2 we show speed tests comparing the 
performance of the MCAMC algorithms to standard MC 
and CT on East model simulations. Fig. 1 shows the tem- 
perature dependence of the CPU time required for gener- 
ating an equilibrium trajectory of total Monte Carlo time 
10 7 x e 2/3 in an East model of N = 10 5 sites. In an s = 1 
CT algorithm the average CPU time for such a simula- 



FIG. 2: System size dependence of CPU time (relative to that 
for CT) for equilibrium East model trajectories of MC time 
( = 3x W 12 /N at T = 0.2. 



tion is independent of T. Fig. 1 shows that at very high 
temperatures standard MC is the fastest method, but as 
T is lowered CT soon outperforms it. At lower temper- 
atures s = 2 MCAMC becomes more efficient than CT 
by a approximately a factor of e 2/3 . As the temperature 
is dropped further, s = 7 MCAMC provides a further 
improvement of approximately e 2 ^ , and so on. 

As discussed above, the efficiency of MCAMC depends 
on the system size. In addition to a reduced time step 
this also determines the probability of encountering the 
isolated entry state for the s > 1 algorithms. Fig. 2 
shows the CPU time, now for different system sizes, at 
fixed temperature and total MC time t = 3 x 10 12 /N. 
Again, the CPU time for such a simulation using CT is 
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FIG. 3: Persistence time r a as a function of inverse tempera- 
ture P = 1/T in the East model in dimensions d = 1 — 5, 10, 13. 
The lines through the data points are quadratic fits, log r a — 
do + aiP + «2/3 2 , with cii fitting parameters. The fit suggests 
that oi ~ 1 in general. The bottom-right panel shows as 
a function of d. This coefficient seems to go as ai ~ b/d, 
with the constant b « 0.8 (shown as a full line). This value is 
between In 2 (dotted line) and (ln2)/2 (dashed line). 



constant. As expected, Fig. 2 shows that as TV becomes 
larger the MCAMC algorithms are less and less effective; 
beyond Ne 2 « O(l) the CT scheme works better. This 
means that in order to maximise the MCAMC efficiency 
one needs to simulate the smallest possible system sizes. 
This is limited by the need to be compatible with bulk 
behaviour, which in the case of facilitated models requires 
that the system in average contains a sufficient number 
of excitations, i.e., Ne cannot be too small. 



IX. EXAMPLE OF RESULTS 

In this section we present an example of numerical re- 
sults obtained with the MCAMC. 

A useful correlation function to study the relaxation 
of facilitated models is the persistence function P(t), e.g. 
ESIHinL which gives the probability that a site has 
not changed its state up to time t. In terms of the lo- 
cal persistence field Pi(t) = 0,1, where 1 indicates that 
site i has not flipped up to that time, and that it 
has flipped at least once, the persistence function reads 
P(t) = N- 1 J2 t Pi(t)- In contrast to standard MC or 
CT simulations, the MCAMC algorithm could run into 
problems when trying to measure persistence. By con- 
struction, it misses some of the events that could occur 
whilst in the transient subspace, for example, from the 
isolated state many spins could flip up and then flip back 
down before finally exiting to an absorbing state. At low 
temperatures in equilibrium, however, the contribution of 
these events is negligible, as the vast majority of changes 
to the global persistence function is from existing exci- 
tations spreading out into unmoved territory, and this is 
captured by the MCAMC algorithm. In fact, the only 
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FIG. 4: Concentration of excitations c(t) as a function of 
scaled time T In i, in the d = 1 East model after a quench from 
infinite temperature, from simulations with s = 7 MCAMC. 



sites one needs to be concerned with are those immedi- 
ately next to the initial excitations, which are very few 
in equilibrium at low T. 

Figure 3 shows the equilibrium persistence time [TR 
Il9j . t q , of the East model in various dimensions d, cal- 
culated using the MCAMC algorithm with s = 2. For all 
dimensions studied we find that r a is a super- Arrhenius 
function of T. This seems to indicate that the East 
model is a fragile in all dimensions. Given that any sim- 
ple mean-field estimate of the relaxation in this model 
would give Arrhenius behaviour, the above result would 
suggest that the East model has no upper critical di- 
mension to its dynamics |l6j| . The data is compatible 
with logT-Q, = ciq + cti/3 + ci2 /3 2 , as expected if relax- 
ation processes in the East model in any d are quasi 
one-dimensional ^(|. The coefficient ai of the quadratic 
fits is compatible with a,2 ~ b/d 16j, with the constant 
b obeying In 2 > b > (ln2)/2, reminiscent of the rigorous 
d = 1 result of Ref . [2C| . Note that the timescales reached 
with the MCAMC in Fig. 3 are between three and five 
orders of magnitude longer than in previous studies |17| . 

The MCAMC proves also useful when simulating out- 
of-equilibrium dynamics. Consider the aging of the East 
model following a quench from infinite temperature. As 
the system relaxes towards its equilibrium the dynam- 
ics proceeds by stages characterized by the distance be- 
tween isolated excitations |2l|. These domains grow as 
d ~ t T ln 2 . Consequently, the isolated transient state also 
plays an important role in such out-of-equilibrium dy- 
namics of the East model, and the MCAMC algorithm is 
also applicable in this regime. Figure 4 shows the aging 
of the concentration of excitations, c(t), with time after 
a quench to low temperatures in the East model, using 
the s = 7 MCAMC algorithm. 

In these aging simulations the nature of the speed-up 
due to the MCAMC becomes evident. Each stage of the 
dynamics is associated with an isolated domain of the 
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form 10 • • • 0. The fc-th stage corresponds to domains of 
typical length I ~ 2 k , and a corresponding energy barrier 
of size k to further relaxation [2l]]. In essence, at each 
successive plateau of c(t) one requires an algorithm that 
produces time steps comparable to the activation time, 
e k @ . An s = 2 MCAMC enables one to push simulations 
one plateau further than CT (s = 1), the s — 7 algorithm 
helps overcome the next energy barrier, and so on. 

X. DISCUSSION 

We have shown that the method of MC with absorb- 
ing Markov chains, MCAMC, of Novotny |(J can be used 
to dramatically speed up simulations of facilitated spin 
models of glasses, such as the East and FA models. Even 
the simplest s — 2 algorithm can improve simulation 
times at low temperature by a factor of e 213 over the n- 
fold or continuous time MC. By increasing the number 
of transient states s even larger computational gains can 
be achieved. One could imagine an algorithm where the 
number of transient states is variable, and at each iter- 
ation the s with the largest number of transient states 
allowed by the current configuration, the smallest dis- 
tance between excitations, is used. 

The next future step would be to adapt MCAMC al- 
gorithms to other interesting KCMs, such as constrained 
lattice gases 0, |H IHEllMl^ and /-spin facilitated 
FA models with / > 1 |7llStl27j . Several features of these 
systems make the application of MCAMC less straight- 



forward: since their kinetic constraints depend on more 
than one site, i.e. facilitation by two or more excita- 
tions in the FA models or two or more vacancies in the 
lattice gases, for generic entry states the tree of possi- 
ble transient states is much larger than for, say, East 
models. This means that the computational cost of the 
necessary bookkeeping will be much higher (bookkeeping 
could be simplified by reducing the possible entry states, 
at the expense of triggering less frequently the MCAMC). 
This problem is compounded by the fact that / > 2 FA 
models are very slow even at moderate temperatures, so 
that the potential exponential in (3 gains from excitation 
rates are very modest, and may not even be enough to 
offset the bookkeeping cost — in constrained lattice gases, 
where barriers are purely entropic, this is bound to be 
worse. In any case, given that the high density or low 
temperature dynamics of these systems is in general so 
much slower than that of East models, a clever MCAMC 
algorithm which overcomes these hurdles could prove ex- 
tremely useful. 
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